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Abstract. 

A hybrid system that combines the advantages of a superconductor flux qubit and 
an electron spin ensemble in diamond is one of the promising devices to realize quantum 
information processing. Exploring the properties of the superconductor diamond 
system is essential for the efficient use of this device. When we perform spectroscopy 
of this system, significant power broadening is observed. However, previous models 
to describe this system are known to be applicable only when the power broadening 
is negligible. Here, we construct a new approach to analyze this system with strong 
driving, and succeed to reproduce the spectrum with the power broadening. Our 
results provide an efficient way to analyze this hybrid system. 
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l. Introduction 

The properties of different types of qubits have been studied extensively over the past 
decade toward realizing quantum computation. Each system has its own advantages 
and disadvantages. The idea of combining different systems is pursued by several 
groups, which would inherit the advantages of each system, such as atoms coupled 
to optical cavities, spins coupled to resonators, and spins coupled to superconducting 
qubits P El El SI El CUE]. 

One of the most successful hybrid schemes is to couple the electron spin ensemble 
with superconducting circuit to realize the fast processing and long storage times. 
Although an electron spin ensemble may possess excellent coherence properties, the 
weak non-linearity makes it difficult to realize a large scale quantum computer. 
Superconducting qubits are sensitive to the noise which leads to a relatively short 
coherence time. Instead, due to a nonlinearity that Josephson Junctions has, a 
superconducting circuit has excellent controllability of quantum states involved in the 
processing, such as single qubit gates, two qubit gates, and projective measurements 
p uni CEH H21 UHl ini HU HZ]. There are two major schemes that utilize this type of 
hybrid system. One of them consists of a spin ensemble, a superconducting resonator, 
and a superconducting transmon qubit |lHl [T9], j2U: ; 2J., [22]. Here, the resonator plays the 
role of a quantum bus [8], [14] where the resonator mediates the interaction between the 
spins and a superconducting qubit. The other one is to use a direct coupling between a 
spin ensemble and the superconducting flux qubit (FQ) [D EH] GEE, 25l l26j. 

A particularly attractive spin ensemble for such hybrid systems is negatively charged 
nitrogen vacancy centers (NV - ) in the diamond crystal whose coherence time is around 
0.6 s [23 EE]. An NV center is a defect consisting of a substitutional nitrogen atom 
and an adjacent vacancy where an additional electron is trapped. It is known that an 
NV - center is a spin-1 system, and there is an energy splitting of 2.88 GHz between the 
states |0) and | ± 1) without an applied magnetic fields. Furthermore, the NV” centers 
can be strongly coupled with superconducting circuits by a super-radiant effect [18, T9], 
where the coupling strength is enhanced by a factor of a /N. 

Several groups have already demonstrated the hybrid coupling systems by the 
ensemble of electron spins and superconducting circuit resonator [201 EH G21 [291 EOl 

m, 32] . Also, Further, Kubo et al. [20] has succeeded in indirectly coupling the 
ensemble to the transmon qubit using a frequency tunable resonator to mediate it. 
However, in order to enhance the coupling strength, this system requires a large number 
of electron spins (~ 10 12 ) in the ensemble to achieve the coupling strength of tens of 
MHz m 2D [221 291 EH SD E2]. These experiments need the mm-size resonator, which 
could be difficult for the integration of many memories on a single quantum based chip 

[2D [22]. 

The advantage to use the spin ensemble and the superconducting flux qubit 
[D 1231 M is that the superconducting flux qubit can couple directly and more strongly 
with the spin ensemble due to the persistent current of the flux qubit around 300 nA 
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~ 900 nA [331 (31] • This system requires a much smaller number of NV - centers spins 
(~ 10 7 ) to realize a strong coupling regime, and the sample occupies only a surface 
area of tens of micro-meter square. Coherent strong coupling between a flux qubit and 
nitrogen-vacancy (NV - ) centers has already been demonstrated [23j [24l 25] , as well as 
the coherent exchange of a single quantum excitation [31] • Quantum memory operations 
involving single qubit and entangled states have been implemented [23]. 

Some theoretical models have been investigated to describe such hybrid systems 
HUIHniM- Traditionally, Jaynes-Cummings model (JC model) is used to describe 
the coupling between the flux qubit and the spin ensemble [7] HB1 HU, 35]. Under a 
strong external magnetic field, the degenerate state |1) and | — 1) of the NV - center 
would be separated far from each other and the NV - center can be considered as a 
spin 1/2 system (two-level system). If the number of excitations in the hybrid system is 
much smaller than the number of NV - centers, we can treat this system as a harmonic 
oscillator and then use JC model to describe the coupling between the flux qubit and 
NV - centers. However, when the applied magnetic field is weak, NV - center will show 
the properties of a spin 1 system (three-level system) [25] because of the degenerate 
states of |1) and | — 1). Therefore, the conventional model cannot be straight-forwardly 
applied to analyze the properties of flux qubit coupled with NV - center for this regime. 

From the spectroscopic measurement in a superconductor diamond system without 
an applied magnetic field, three resonant peaks have been observed in [2T[ [31]. A sharp 
peak has been clearly observed in the middle of the avoided crossing caused by the 
coherent coupling between the superconducting circuits and an ensemble of electron 
spins in nitrogen-vacancy center. The center narrow peak shows a longer lifetime than 
that of the other two broader peaks, which may provide an alternative approach for the 
quantum memory [25] [36]. 

In the paper [25], a full Hamiltonian model for the hybrid system of the flux qubit 
and NV - centers is used to interpret the mechanism causing the sharp narrow peak, 
which contains the properties of spin 1 system such as inhomogeneous strain, randomized 
magnetic field from PI center (substitutional nitrogen centers in diamond), and zero field 
splitting fluctuation [25] . The NV - center is regarded as a three-level system, and it was 
shown that two collective modes of the NV - ensembles are relevant for the spectrum. 
The two broader peaks are associated with the dressed bright state in NV - that can 
be directly coupled with the flux qubit. The center narrow peak arises from another 
specific collective mode that cannot be directly driven by the flux qubit. This is called 
a ’’dark state”. Due to the effects of the strain and Zeeman splitting, the bright state 
and dark state can exchange the excitation, while the bright states interact with the 
flux qubit. In that theoretical analysis [25], the flux qubit is regarded as a harmonic 
oscillator in the limit of low excitation energy so that the Heisenberg equations can be 
solved in the frequency domain. However, this approximation is valid only when the 
driving power is weak, and will be invalid to describe the phenomenon when the power 
broadening is observed in the strong driving power case. Since it is not always possible 
to obtain a reliable data without power broadening due to a small signal to noise ratio, 
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this could be a limitation to describe the hybrid system for such a regime. Therefore, 
it is necessary to build an alternative approach that can include the effect of power 
broadening in the strong driving regime. 

In this paper, we extend the previous approach to reproduce the experimental 
spectroscopy that is significantly affected by the power broadening. The rest of this 
paper is organized as follows. In section 2, we introduce the previous model for this 
hybrid system [25]. In section 3, we highlight a new theoretical model to show that the 
numerical results agree with the experiments. Finally, section 4 contains a summary of 
our results. 


2. Hamiltonian 


Let us review the Hamiltonian introduced by the paper in [25]. The diamond containing 
NV” centers is attached on top of the flux qubit so that these two system can interact 
each other, as described in the Appendix 15.11 The Hamiltonian of this system includes 
not only the interaction between the flux qubit and NV” centers, but also the strain 
distributions, the inhomogeneous zero held splittings, the effect of PI centers, and the 
hyperhne coupling from the nuclear spin, described as follows: 

P //flux T Pdrive T Pens T Pint (1) 


The flux qubit Hamiltonian, interaction Hamiltonian, the NV 
Hamiltonian and the driving Hamiltonian can be written as: 

TT h „ h 

Pfiux T —/\o'x 

Pdrive = flX COS Ult ■ CT Z 
N 

if. n, = n V { D ^lk + Wh - sit) + 

k =1 

Ei k) (S x , k S yjk + S y>k S Xtk ) + MbBS • s k } 

N 

Pint = hge^BCTz ( ^ B^ ‘ Sfc) 

k =1 


diamond ensemble 


( 2 ) 

(3) 

(4) 

(5) 


where a XjZ denotes the Pauli matrices for the flux qubit, a z represents the population 
difference between two persistent current states in the flux qubit, A denotes the flux 
qubits tunneling energy, and £ denotes the bias. From the Hamiltonian of the flux 
qubit, the frequency of the gap tunable flux qubit is calculated as ujfq = Ve 2 + A 2 . 
Since the bias e can be controlled by the external magnetic flux, we use the case of 
no bias (where the bias e is equal to 0), in which case the flux qubit has the best 
coherence time at the optimal point [22, 24] . In the expression for the driving field of 
the flux qubit P drive, X is the amplitude of the field, and c o is the angular frequency of 
the microwave. H ens represents the electron spin ensemble of N individual NV” centers 
with S Xt y tZ describing the electron spin 1 operators of the individual NV” center. H ens is 
characterized by the zero-field spitting P, the strain induced spitting E and the Zeeman 
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splitting term g e g, B B N v ■ S. Here, g e = 2 is the NV” Lande factor and fi B is Bohr’s 
magneton. B^v is the magnetic field composed by the inhomogeneous magnetic field 
from PI centers and hyperhne induced magnetic held generated by nitrogen nuclear 
spins. E[ k> denotes the strain along the x direction while denote the strain along 
the y direction. The term B q b in the Hi nt presents the magnetic held generated by a 
hux-qubit persistent current. 

Since it is difficult to solve the Hamiltonian described in the Eq. [t] our previous 
model considered the hux qubit and NV” centers as harmonic oscillators on the condition 
that the driving held strength of the hux qubit is weak [25]. We call this a many 
harmonic oscillator model (MHOM). The details of this model are discussed in the 
Appendix 15.31 Moreover, to obtain an intuitive explanation for the spectroscopic 
result, the Hamiltonian is further simplified with homogeneous NV” centers under a 
zero magnetic held (such as D *. = D, B^ = B and = E ) without the driving term 
(A = 0) [25]. Considering a homogeneous system, one can describe the NV ensemble by 
using just two harmonic oscillators. The Hamiltonian is rewritten as: 

H = huFQC^c + Hujb&b + hu>d<ftd + hg{$b + cb t) + hJe l6 b^d + hJe~ l6 b$ ( 6 ) 

where u>b (cod) denotes the frequency of the bright (dark) mode of the NV” centers, $ 
denotes a creation operator of the hux qubit, ¥ (cr) denotes a creational operator of 
the bright (dark) mode of the NV” centers respectively, g denotes the coupling strength 
between the NV” centers and the hux qubit, Je ld = g e ^BB + %E denotes the coupling 
between the bright mode and dark mode. It is worth mentioning that the bright mode 
can be coupled with the hux qubit while the dark mode has no direct coupling with the 
hux qubit. The frequency of the bright state is assumed to be equal to that of the dark 
state, and so we have uj^v = u>b = = D. This corresponds to the zero-held splitting 

of the NV” center between |0) and | ±1) states at zero magnetic held. (In the Appendix 
15.41 the eigenvalues and eigenvectors of this simplihed Hamiltonian are described, and 
we discuss the properties of them.) 

3. Main results 

3.1. Master equation model (ME) 

A spectroscopy of this system was performed, and it was shown that power broadening 
due to the strong driving power is relevant in this system (shown in the [25]). There 
are several models attempting to understand this system. However, no existing model 
can explain the power broadening induced by strong microwave pulses for this system 
[ 251 [3-4] . It is important to build a model to understand this phenomenon. Such a 
model would contribute to evaluate the performance of the hybrid system under different 
driving powers and would help to optimize the fabrication parameters. 

Here, we introduce our approach using a master equation model (ME) to study 
this hybrid system (shown in Figfl]). As the ensemble is driven by the flux qubit that 
is an effective single photon source, the number of excitations in the ensemble is much 
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smaller than the number of electronic spins, and therefore it is reasonable to regard the 
spin ensemble as a number of harmonic oscillators as the previous authors considered. 
However, when the external driving power for the flux qubit is strong and induces the 
effect of the power broadening, we should consider the flux qubit as a two-level system 
and use the spin 1/2 operator a instead of the harmonic oscillator operator. This is 
because strong driving power would populate higher energy structure in a harmonic 
oscillator system, which is different from the two-level system. 


Hybrid system 



Figure 1. Schematic of the model for our approach (ME) to describe superconductor 
diamond hybrid system. The flux qubit is considered as a two level system. On the 
other hand, the NV centers are replaced by two harmonic oscillators, which we call a 
bright state mode and a dark state mode. The bright state is coupled with both the 
flux qubit and the dark state. 

In such a case, the Hamiltonian should be given without harmonic oscillator 
approximation: 

H = ~ u)d z + h(coNV ~ oj)b^b + h(ujNv ~ co)$d + hg(a + b + 

+ hJe ie b ] d + hJe~ i0 bd f + ^a x (7) 

Moreover, we use the dissipative term of Lindblad equation [371 ESI [39] to present the 
effect of inhomogeneous broadening induced by the strain and Zeeman splitting from 
each NV~ center: 

! = -i[p,H] + L(p) (8) 

L(p) = - Y FQ (a+a_p + pd+oL - 2dLpd+) 

— Y b (b^bp + p&b — 2 bpP) — Y d (Sdp + p$d — 2 dp$) (9) 

where L(p) denotes a Lindblad super-operator and p is the density matrix. Y F q, T ft , and 
r d denote the decay rate of the flux qubit, bright state, and dark state, respectively. The 
inhomogeneous broadening induces unknown phases on the bright state (or dark state), 
which transforms the collective mode into another orthogonal modes [24^ :25l [34] . It has 
been shown that an inhomogeneous distribution of the frequency can be considered as 
an energy relaxation of the collective mode [35], which can be described as a Markovian 
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process. Actually, an exponential decay of the collective mode in this hybrid system has 
been observed in [23], which is consistent with the model of the Markovian noise. So, 
to include the effect of inhomogeneous broadening, we add the Markovian relaxation 
terms of the bright state and the dark state in the Lindblad equation. 


3.2. Three Harmonic Oscillator model (THOM) 


In order to describe our sample of the flux qubit coupled with the NV” centers, we 
have already obtained the necessary parameters (£, SB, 5E, SD) for the simulation of 
the MHOM from the experiment [253. However, as we describe before, this model will 
be invalid for a strong driving regime, and we should use the ME when the power 
broadening is relevant. Since we use a different model, we cannot straightforwardly 
use the parameters that are used in MHOM. We need to find suitable parameter 
set (g, J, Tfq, T*,, T^) for the ME. For this purpose, it is convenient to have a fitting 
function that let us End out such a new parameter set (g, J, T F q, T b , T d ) from the 
known parameter set (£,5B,5E,5D). In the condition of weak excitation, the excited 
probability of the qubit is described as 


(a+a.) 


V 2 ,_ [(iU-HKirj-M')- J 2 ] _ 

2 ’ 1 (ir c - 4)[(;r t - J t )(iT d - w ') - ,P] - gH,r d - w ') 


where uj' F q = lu F q — lu and lu' nv = ~ w. We can derive this from a three harmonic 

oscillator model (THOM) where the flux qubit is regarded as a harmonic oscillator (See 
the Appendix 15.51 for the details). It is worth mentioning that this analytical solution 
is useful only for weak driving limit. However, this approach plays a crucial role to find 
the parameter set for the master equation approach, as we describe later. 


3.3. Determination of parameters for ME 

In this section, we will introduce the process about how to determine the necessary 
parameters to describe our devices for our ME. We have already obtained the parameters 
to describe our devices for the MHOM (the model described in the paper [25]): 

SB/2tt = 0.056 rnT; <5E/27 t = 4.4 MHz; 5D/27T = 0.2 MHz; (11) 

We construct a scheme to translate the parameters (£, SB, SE , SD) in the MHOM into 
the parameters in our new approach (g, J,T b ,r d ), which is one of the main results in 
our paper. As mentioned before, THOM provides the fitting function to obtain such 
parameters. If we fit the MHOM spectrum by this function, it seems that we might 
determine the parameters needed in the ME because THOM is equivalent to ME for a 
weak driving regime. However, due to many parameters to be estimated, if we naively 
apply the fitting function to the MHOM, we cannot obtain a unique set of parameters 
from the fitting result. In general, it requires other methods to fix the set of parameters, 
and we have developed such method. Firstly, from the T 1 measurement, we estimate 
the value of the flux qubit decay rate T F q. It is worth mentioning that T F q denotes a 
decay rate of the flux qubit for both MHOM and ME so that we could use the same 





value for T FQ in the ME as that used in MHOM (T fq — -^r — 0.33 MHz). Secondly, 
by using the expression of eigenvalues and eigenvectors of the Hamiltonian described in 
Eq. [?1 we End the method to estimate the values of J and g from MHOM. Finally, by 
fixing the values of J, g, and T F q, we use the analytical solution of THOM to fit the 
MHOM’s result and determine the other two parameters of W and T^. The detail of 
the procedure is described as follows. 

3.3.1. Estimation of the coupling strength J and g We explain how to determine the 
coupling strength J and g from the MHOM. 

First, by the Hamiltonian in the Eq. [7J the frequency difference between the right 
peak and left peak on the resonant condition depends on the coupling strength J and g 
( 5 = 0 ), as shown in the Appendix 15.41 The relationship is described as 

{Eright — Ei ef t)/h = 2 a y g 2 + J 2 (12) 

Also, if we assume g T^,, T^, T c in the Eq. OH we can calculate the frequency difference 
between the left and right peak 

(E right - E left )/h ~ y/ 4(g 2 + J 2 ) - (T b - T d ) 2 ~ 2\Jg 2 + J 2 (13) 

and this is consistent with the Eq. [HI Moreover, from the frequency difference between 
the left and right peak plotted by MHOM (shown in Fig. [8]), we obtain the following 
relationship 

(■ E right - Eieft)/h ~ 27 MHz A/g 2 + J 2 ~ 13.5 x 2vr MHz (14) 



Figure 2. Energy shift of the middle peak in the spectroscopy against the detuning. 
Here, the red dots represent the simulation result with the MHOM, and the green dots 
represent numerically calculated eigenvalues. The blue line represents an analytical 
solution using a perturbation theory. 

Second, we use a relationship between the middle peak frequency and detuning 
6. We have shown that the middle peak frequency of the spectroscopy will be shifted 
by adding the energy detuning, and the energy shift increases linearly as we enlarge 
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the detuning when the detuning is considered as a perturbation (5 <C ujnv) (See the 
Appendix 15.51) . This is described as 


SJ 2 


^middle ^NV — 


g 2 + J 2 


(15) 


where the proof is described in the Appendix 15.51 Actually, we plot the energy shift of 
the middle peak against the detuning by solving MHOM (shown in FigJSJ), and we fit 
this plot by Eq. [15j From this fitting, we obtain the following relationship between 
the parameters for ME: 

j 2 

—-- « 0.067 => J ~ 0.27 g (16) 

g 2 + J 2 


However, since we use a perturbation to obtain the fitting function in EqJl5l this result 
would be invalid for a large detuning. To confirm the validity of the perturbation, we 
numerically calculate the eigenvalues of the Hamiltonian in the Eq. [71 and plot the 
results in Fig. [2} Since we have a good agreement between perturbation calculation and 
numerical results, we conclude that the perturbation is valid in this parameter regime. 

By combining the Eq. [14] and Eq. [T6l we can estimate the values of J and g: 


J ~ 3.5 x 2-7T MHz; g ~ 13 x 27T MHz; 


(V) 


3.3.2. Estimation of decay rate T^ and The value of T^ mainly affects the width of 
the side peaks while the width of the middle peak is determined by the value of (See 
the Appendix 15.51 for the details). This means that, we can use the THOM as a fitting 
function to fit the spectrum reproduced by the MHOM if the parameters T^q, J and g 
are known, and we have obtained the following parameters 

~ 0.49 x 27 t MHz; Tb — 6.4 x 27 t MHz; (18) 

Therefore, all necessary parameters for the ME have been estimated by our scheme. 


3 . 4 . Reproducing experimental results with the ME 


With these parameters, we use the ME to reproduce the spectroscopic measurement 
in experiment for different driving power (See Fig[3]and Fig[4j). Not only in the weak 
power case but also in the strong driving power case, there is a good agreement between 
simulation and experiment, as shown in Fig[3] and FigJTJ Thus, our new approach 
with the ME is shown to be useful to reproduce the spectroscopy even when the power 
broadening is relevant. 

Furthermore, in order to compare the theory and experiment in detail, we analyze 
the width of the middle peak which can be well fitted by a Lorentz function. 

„2 

+ c (19) 


ay 


(oj — a /) 2 + y 2 

where 7 correspond to the width of the Lorentz curve and uj' is the center frequency, a 
and c and the fitting coefficients. We use the Lorentzian function to fit the middle peak 
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a) Driving powei^-16dBm 



D riving frequency (G Hz) 



Driving frequency (GHz) 


Driving power=-20dBm 




Figure 3. Spectroscopy with different driving power. The blue dots are the measured 
data where continuous lines are drawn through the points as a guide to the eye. The 
brown line denotes the error bar. The red curve shows the numerical results by solving 
the ME. For the simulation, we use parameters of lufq = <^>nv — 2.878 x 2-7T GHz, 
g = 12.95 x 2-7T MHz, J = 3.46 x 27r MHz, T fq = 0.300 x 2-7T MHz, W = 0.493 x 2-7T 
MHz, and T{, = 6.433 x 27r MHz. Here, a, b 7 c, and d present the spectrum in different 
driving power: —16 dBm, —20 dBm, —24 dBm and —30 dBm, respectively. It is worth 
mentioning that we cannot observe a power broadening around -30dBm as described 
in Fig. [5j and so the power is sufficiently weak in this regime. 


calculated by the ME (shown in Figj5]a). By the fitting, we obtain the FWHM (Full 
Width at Half Maximum) 27 with the ME, MHOM, and the experiments, respectively. 
In the Fig j5j)b) , we plot them, and the power broadening is clearly observed in both 
the experiment and the ME, while MHOM cannot describe the power broadening. 
These results also show that the ME surpasses the MHOM in order to reproduce the 
experiments. 

4. Conclusion 

In summary, we introduce a new theoretical approach to describe spectroscopic 
measurements of a superconducting flux qubit coupled with NV - centers in diamond. 
Although previous models are applicable only when the driving power of the applied 
microwave is weak, we have succeeded to reproduce the experimental spectroscopy even 
when the power broadening becomes relevant due to the strong driving. Since it is 
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Figure 4. The spectrum with power broadening, (a) Power dependence of the energy 
spectrum of our hybrid system obtained from the experiment, (b) Numerical results of 
the energy spectrum by solving ME. (c) Numerical results of the energy spectrum by 
solving MHOM. We use the same parameters as those in FigEl Inset of each picture 
is the enlargement of the spectrum over the middle peak region. 


a) b ) 



Figure 5. Behavior of the width of the middle peak with different driving power, (a) 
The middle peak is fit by the Lorentz function. The red dots denote the simulation 
result from ME. The blue curve denotes the fitting result of a Lorentz function. The 
parameters for ME are the same as Fig 16. (b) The full width at half maximum 

(FWHM) with different driving power. The blue dots present the FWHM from the 
experiment with error bar. The red line denotes the FWHM from the ME by Lorentz 
fitting. The green line is the FWHM resulted from MHOM simulation. 
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typically difficult to remove the effect of the power broadening in the spectroscopy when 
a superconducting flux qubit is driven by the microwave, our results provide an efficient 
way to analyze the superconductor diamond hybrid system. Our method will be useful 
to characterize this system for the application of quantum information processing. 
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5. Appendix 

5.1. Experimental setup 

We briefly describe the experimental setup. The hybrid system is shown in Figj6] The 
NV - ensemble is generated by an ion implantation and annealing in vacuum |34j. The 
density of the NV - centers is approximately 5 x 10 17 cm -3 . The diamond including 
the NV - centers is glued on the top of the flux qubit where the distance between the 
flux qubit and the surface of the NV - diamond is less than 1 pm. We can control the 
hybrid system by the microwave, and this system is measured by switching current in 
the SQUID that is inductively coupled with the flux qubit. 


> 

' f 

< - 


Control line for the FQ 
gap tuning 



< - 

' f 

-- 

Microwave line for 
the driving 


Figure 6. Illustration of the hybrid system. It contains the SQUID, the control line, 
microwave line, a flux qubit with four junctions (the cross in the diagram), and an 
ensemble of NV center in the diamond that is glued on the top of the flux qubit [34]. 
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5.2. Many Harmonic Oscillator model (MHOM) 

There are two groups that performed spectroscopy of a superconductor-diamond hybrid 
system under zero external magnetic field j2T, 23 04], and both groups observed two 
broader side peaks and a middle sharp peak in the spectrum. However, in the papers 
[21101], such three peaks could not be reproduced by the theoretical model. In a paper 
[25] , such three peaks are firstly reproduced by a numerical simulation based on a new 
model. 


NV Centers 


Bright state 

|+ 1 >+|- 1 > 


Dark state 

|+ 1 >- |- 1 > 



di 

d 2 


-*■ d 


Figure 7. The schematic of the MHOM. Here, c denotes the flux qubit, u>( k \ (w^d) 
denotes the frequency of the fcth bright state (a dark state) in NV - centers, b\. (d\) 
denotes the bright state (dark state) of the fcth NV - centers, J k denotes the effect from 
Zeeman splitting, J k denotes the effect from the strain, and denotes the coupling 
between the flux qubit and fcth NV - centers. 


We explain how to obtain the model introduced in [25]. In order to solve the 
full Hamiltonian described in Eq. [H the flux qubit and NV” centers are regarded as 
harmonic oscillators on the condition that the driving held strength of the flux qubit is 
weak [25]J . In this case, the average number of the excitation at the flux qubit and NV” 
center is much less than one. Since there are two types of excited states in the NV” 
center ensemble, two creation operators (such as b + and d + ) are defined to describe 
the properties of NV” centers (shown in Fig]?]). Here, b + denotes a creation operator 
of the bright states that can be directly coupled with the flux qubit, while d + denotes a 
creation operator of the dark states that interact only indirectly with the flux qubit via 
the bright states. Then making the rotating wave approximation, we define the many 
harmonic oscillator model (MHOM) and obtain the Hamiltonian as following: 

H = h{oj F Q - u)c^c + h^(c + c t ) 

N 

+ ^ [H UJ b k) ~ ~ u)d\d k + h( k (c ] b k + cb\) 

k =1 

“I - fiiyJk iJ k )bkd k 


( 20 ) 
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where u/ k \ = D k — E^ denotes a frequency of the fcth bright state and uj^d = D k +E^ 
denotes a frequency of the kth dark state. Since the effect of the the Zeeman splitting 
(strain) allows the transition between the bright states and dark states, a coupling 
strength of J k = g e fi B B k (j' k = E^) is defined to present the interaction between them. 

Similarly, the two-level system of the flux qubit is simplified as a harmonic oscillator, 
and so this can be described by the creational operator c + . The flux qubit is coupled 
only with one mode of the ensemble (the bright states) and the individual coupling 
strength is presented by the value of (. A denotes the driving power of the flux qubit 
while u denotes the frequency of the driving power. 


5.3. Heisenberg-Langevin equations of MHOM 


We explain the details about how we can solve the model introduced in [25], which we 
call MHOM. Based on the Hamiltonian in Eq. [20j the Heisenberg-Langevin equations 
of the system are written as |25]: 


£-c = -{Y FQ + iu FQ )c-i(^Ck- b^j 

k =1 

fb k = - i C, I - (ij k + j' k )d k -i( k -c 

fd k = - (r d + - (iJk + J' k )b k 


( 21 ) 

( 22 ) 

(23) 


where T F q, T t,, denote the energy decay of the flux qubit and the NV - center bright 
state and dark state, respectively. By transforming these into Fourier space, then we 
obtained the expression of c: 


^ A 
C “ 2 X 


w - ufq + iV FQ - ^ Y^k=i Kfcl 2 




+*rd 




(24) 


These Heisenberg-Langevin equations was solved to obtain (c^c), which denotes 
the average number of the excitation in the flux qubit [25]. It is worth mentioning 
that the main decoherence source of this hybrid system is inhomogeneous broadening 
[211 [221 E3l [2H ]20l |29l [30l EU |32], and the Heisenberg-Langevin equations used in [25] 
include such effect. In this model, the spectroscopic measurements of the flux qubit 
with NV - centers are reproduced with a narrow peak located in the middle of the avoid 
crossing (shown in Fig. [8]). Two broad resonances are observed on the side, which we 
call a left peak and a right peak, respectively. Also, one sharp resonance is observed 
between them around 2.88 GHz, which we call a middle peak. These agree with the 
experimental results in the limit of weak driving [2T[ 2.5. 154] . 


5.f. Eigenvector and Eigenvalue 

To understand the spectrum shown in Fig. [8] the Hamiltonian in Eq. [6] is diagonalised 
in [25], and we explain the results here. 

H = hu> F QC^c+hujNv&b+fruj F rvd ji d+hg(c l b+c& i )+hJ(e ie b^d+e~ ie b$) (25) 
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Figure 8. The energy spectrum of the hybrid system by MHOM in [25j. Here, 
x axis denotes the microwave driving frequency and y axis denotes the switching 
probability of the SQUID which corresponds to the microwave absorption of the flux 
qubit (Pswitching = 1 — Pabsorption)■ The spectrum is plotted under the conditions 
of g e fj,BB/2ir = 28 MHz; N = 36000; u}fq/2it = Ub/2n = ujd/2n = 2.88 GHz; 
S(g e g B B/2 tt) = 3.1 MHz(FWHM); 5(E/2n) = 4.4 MHz(FWHM); T fq /2tt = 0.3 
MHz(FWHM); T NV /2n = S(D/2 tt) = 0.2 MHz(FWHM); X/2n = 20 MHz. 


where g denotes the collective coupling strength between the bright state and the flux 
qubit, J denotes the coupling strength between the bright state and the dark state, and 


FQ |t> 


Oi F Q 


IQ 


FQ coupled with N V 



Figure 9. The energy diagram. In the picture, the left side denotes the energy 
diagram without coupling (coupling strength g = 0) while the right side denotes the 
case with coupling. Here we assume that the coupling strength between the dark state 
and bright state is much smaller than the coupling strength between the flux qubit 
and NV - center (J <C g). 


|0) denotes the ground state of NV center. Since NV center is a spin 1 system, it has 
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three electron states: |0), |1) and | — 1). Suppose that | B) and | D) denote the excited 
states of bright state and dark state respectively where 


|B) = -h(|l) + |-l» 
|D> = 4(|1)-I-1». 


(26) 

(27) 


| J,) and | t) denote the ground state and excited state of the flux qubit. The relationship 
between the ground state and the excited state is described as follows: 

^|0)w = \B) nv (28) 

<^|0)w = \D)nv (29) 

C f | 1)fq = I T)fq; (30) 


The transition from the ground state to the first excited state (shown in Figj9]), which 
corresponds to the left side peak in the experimental spectroscopy (shown in FigE]), is 
described by: 


A Eleft = Eleft — E 0 = h(u>NV ~ VQ 2 + J 2 ) 


I Euft) = ( 


7 I 1 Bi) + 


1 9 

V2^TT* 


lot)) + 


1 Te~ ie 

-p . ior) 
y/2 


(31) 

(32) 


The transition from the ground state to the second excited state (shown in FigEJ), which 
corresponds to the middle peak of the spectroscopy (shown in FiglH]), is described by: 


A E rnu ldlc E m iddl e Eq h/jJ J\y 

Je ie 


J middle 


) = 


=|ot) + 


v / 7+4 i 


I Di) 


(33) 

(34) 


The transition from the ground state to the third excited state (shown in FiglH]) , which 
corresponds to the right side peak of the experimental spectroscopy (shown in FigJH]), 
is described by: 

A E r ighf Eright Eq ll(uj]yy -\- T J ) (35) 

'*"•»“) = ( 7! |b w + 71 1}) + 717& |d 4> (36) 


Due to the hybridization induced by the coupling between the flux qubit and NV - 
centers, the energy level of the exited states of the hybrid system is split into three 
levels. The left and right peaks in the Fig.((8]) correspond to the first and third excited 
states which contain a bright state of the NV - centers. We call these two exited states 
as the hybrid bright states. The middle peak correspond to the second excited state 
that shows the existence of a NV - dark state, which we call the hybrid dark state. 
Interestingly, the second excited state contains |0 f) where the flux qubit is excited so 
that the signal of this excited state can be detected via the spectroscopic measurement 
of the flux qubit. 
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5.5. Properties of THOM 

We explain how to solve THOM, and discuss the properties of this model. 

5.5.1. Solution of THOM In the condition of weak excitation, we can use a harmonic 
oscillator to replace the spin 1/2 operator (<r + cf) and use two harmonic oscillators O 
and df to represent NV~ center. Then, we use the non-hermitian decay terms ( iT c Oc 1 
iTiffb, iTddfd ) to take the role of Lindblad terms, presenting the relaxation in the flux 
qubit, bright state and dark state (shown in Fig JTOl). We rewrite the Hamiltonian as 


Hybrid system 



A(co) 


NV center ensemble 

a } /Tn 


Microwave Flux qyibit Bright state Darkgtate 
driving power V _A__N_ Ay 


1 FQ 


Figure 10. Schematic of the model for the THOM to describe superconductor 
diamond hybrid system. The flux qubit is considered as a harmonic operator. The 
bright state is coupled with both the flux qubit and the dark state. 


follows. 


H = hjj F q(fc + huiNvb'b + hujjyydf d + hgfcfb + cH) 

+hJ(e l6 b^ d + e~ ie bd^) + A(c/ + c) cos (cot) — iT c (fc — iT b lfb — iT d d)d (37) 


, which we call a three harmonic oscillator model (THOM). Here, A denotes the 
amplitude of the driving power and oo denotes the frequency of the driving power. T 
denotes the energy relaxation term for flux qubit, bright state and dark state. 

After performing the rotating wave approximation with the Hamiltonian described 
in Eq. (1371) . we have: 


H = Huj F qC^c + hu NV b^b + huj NV $d + hg(c l b + elf) + hJ{e l6 b^d + e l9 b$) 

+-A(c/ + c) — iTfcfc — iT b Pb — iT d dfd (38) 

where uj' F q = lo F q — uj denotes the difference between the flux qubit frequency and 
driving power frequency. u' NV = ~ w denotes the frequency difference between 

NV - center and microwave driving power. 

Based on this Hamiltonian, we obtain the following Heisenberg equation: 



— i[c, H] = —i{ijj F QC + gb+ A — iT c c) 


(39) 
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Figure 11. The relationship between the spectrum and energy level. In the picture 
(a), we plot the spectrum using the THOM where ojfq/ 27t = lonv/^k = 2.88 GHz, 
g = 13.0 x 2-7T MHz, J = 3.46 x 2-7T MHz, T ^q = 0.30 x 2-7T MHz, T^ = 0.50 x 2-7T 
MHz, T;, = 6.40 x 27 t MHz, and A = 1.00 MHz. The spectrum denotes the transition 
from a ground state to each excited state. In the picture (b), we describe the energy 
levels diagram. The transition between the ground state and the first (third) excited 
state | Ei e ft) ( \E right )) corresponds to the left (right) broad peak in the spectrum. The 
transition between the ground and the second excited state \E m i d die ) corresponds to 
the middle narrow peak. 


—b = — i[b, H] = —i(uj NV d + gc+ Je l9 d — iY b b) 

LLL 

—d = — i[d, H] = —i(pj NV d + Je~ l6 d — iT^d) 

LLL 


(40) 

(41) 


We assume that, after driving the system for a long time, we obtain a steady state 
and so we obtain: 


(4«)«oo = = (4<0<-.°o = 0 

dt dt dt 

By solving these, the excited probability of the qubit is described as: 


(0 , c) = (l) 2 x 


[(;r„ - «l)(<r, - u d ) - j 2 ] 


(42) 


! (43) 


(ir c - w;)[(ir b - u' b ){iT d - u' d ) - J 2 } - g 2 {iT d - J d ) 

where J F q = Ufq — OJ and u' NV = ~ u. We plot this in Fig. [1TJ and this can 

reproduce three peaks that are observed in this hybrid system |25j . 

From the THOM, we confirm that the phase 6 in the expression of coupling between 
the dark state and the bright state would not affect the spectroscopic measurement 
results, and therefore, we set 6 = 0 throughout this paper. 


5.5.2. Properties of the THOM with a resonant condition We investigate the behavior 
of this THOM for several parameters when the flux qubit is resonant with the NV - 
centers. From the simulation result, the hybrid system contains three excited states, 
and the middle peak in the spectroscopy denotes the hybrid dark state whose width is 
much narrower than that of two side peaks. In order to understand the effect of each 
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2.84 2.86 2.88 2.90 2.92 

Driving frequency (GHz) 


Figure 12. Density plot of the spectroscopy of the THOM to describe the dependency 
on J and g where T pq = 0.30 x 2-7T MHz, = 0.50 x 27 t MHz, lb = 6.40 x 2-7T MHz, 
and A = 1 x 27 t MHz. In the picture (a), as the value of g increases, the middle 
peak becomes darker. Meanwhile, the two side peaks will be largely split for a larger 
g. On the other hand, in the picture (b), as the value of J increases, the effective 
coupling between the dark state and the flux qubit becomes stronger where the bright 
state mediates the interaction, which makes the middle peak higher than the other two 
peaks. 


parameter, we draw a density plot and show the relationship between the parameters 
and spectroscopy (shown in FigJl2l and FigJT3l). 

As the value of g increases, the separation between two peaks becomes larger (shown 
in Fig{12]a). This can be explained by an increase of the vacuum Rabi splitting that is 
determined by the J and g as follows. 

bright — AEi e f t = 2 h\J g 2 + J 2 (44) 

where the value of g depends on both a persistent current of the flux qubit and the 
distance between NV~ center and flux qubit. The value of J is determined by the 
inhomogeneous effect of NV - centers. Since g is usually much larger than J, we have 
AE right — AEi e f t ~ 2h(g + A'j (here, g J) and so this energy gap is mainly determined 
by the value of g. The area of the peak is determined by the form of the eigenvectors. 
From the expression of \E midd i e ), as g increases, the weight of \D ),) increases while the 
weight of the flux qubit |0 t) decreases. This means that, for a large g, it becomes 
more difficult to detect the dark state by the spectroscopy where only the flux qubit 
excitation is measured, and the middle peak disappears (See the Fig [12] a). 

As the value of J increases, the coupling between the dark state and bright state 
becomes stronger, which induces a stronger hybridization between dark state and flux 
qubit. Since the weight of |0 t) increases with enlarging J, the area of the middle peak 
becomes larger (See the Fig [12] b). 

If we change the decay rate T& of the bright state, it will make two sides peaks 
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a) 


c) 



2.84 2.86 2.88 2.90 2.92 

Driving frequency (GHz) 



D riving frequency (GHz) 


Figure 13. Density plot of the spectroscopy to describe the dependency on the decay 
rate T<j, and Tfq, respectively. We fix the other parameters g = 13.0 x 2i r MHz; 
J = 3.46 x 2tt MHz; A = 1 x 2tt MHz and then change the Ft,, T<j, and T F q values. 
These three pictures (a), (b) and (c) correspond to the excited population varieties 
with T;,, Td, and T fq- In the picture (a), Fb mainly affects the width of the side peak 
which becomes broader as the value of IV increases. On the other hand, the width of 
the middle peak is insensitive against IV In the picture (b), the relaxation rate of the 
dark sate IV has a significant contribution to the excited population of middle peak 
but has a smaller influence on the side peaks. The picture (c) shows the effect of T fq, 
and it change the width of the side peaks and has insignificant effect on the middle 
peak. The side peaks decay rapidly as T fq increase. 


broader (shown in Fig{13]a), due to a shorter life time. For the effect of T^, it will only 
contribute to the middle peak and has almost no influence on the side peaks (shown in 
Fig{j3]b), which is consistent with the form of the eigenvector described in Eql32l and 
Eq]36](J <C g). The flux qubit is directly coupled with bright state and is coupled with 
the dark state indirectly. Therefore, the flux qubit decay rate F F q mainly affects the 
width of the side peaks, while it has smaller effect on the middle peak (shown in FigJl3l 
c). 

5.5.3. Properties of the THOM with a detuning Next we consider a case to add 
the detuning between the flux qubit and NV _ center. We define the detuning as 
S = ojfq — ^nv- As long as the detuning is much smaller than other parameters, 
it is valid to regard the detuning as the perturbation to the system. We rewrite the 
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Figure 14. The energy level diagram with and without detuning. Picture (a) shows 
the energy level and energy gap between ground state and each excited state without 
detuning. Picture (b) shows the effect of inducing the detuning between flux qubit and 
NV - center. From the energy level diagram, the frequency of each peak will cause a 
blue shift depending on the 8 value. 


Hamiltonian as following: 

H = h((jj N y + 5)$c+f\WNv^b+fa^NV$d+hg($b+cV i ) + hj(iJd+b(ft)( / \:b) 

By diagonalising the Hamiltonian and using the perturbation theory, we can 
calculate the eigenvalues and eigenvectors (shown in FigfTdl). 


I g 2 

AEleft = Eleft - E 0 « h(uj NV - vV + J 2 + ^ g 2 + J2 \ 


(46) 


I Eleft) — -7^[1 + 

1 9 


gH 
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I B\) 
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J 2 S 
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1 J r _ g 2 5 J 2 5 


;[ 1 - 


4- 
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AEjniddle -SmieMZe Eq ~ bl(ujNV 4“ & 
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Figure 15. Energy shift of the middle peak against detuning for MHOM. Each picture 
presents for the case of the detuning as 0, 2 MHz, 6 MHz and 10 MHz, respectively. 
The detuning makes it difficult to excite the left and middle peak whose frequency is 
far from the energy of the flux qubit. Here we choose the same parameters as in Fig 
[8] except of the frequency of the flux qubit which is equal to the sum of the detuning 
and frequency of NV - center. 


a) b) 



‘ ‘ ■ 1 ‘ ‘ ■ 1 ■ ■ ■ ■ 1 ■ ■ ■ ■ • ‘ . .. 

2.880 2.885 2 890 2.895 2.900 2.905 2.910 2.880 2.885 2.890 2.895 2.900 


The fluxqubit frequency(GHz) The fluxqubit frequency(GHz) 

Figure 16. Numerical simulation by solving THOM where the x axis denotes the 
flux-qubit frequency and y axis denotes the microwave frequency, (a) Density plot of 
the spectroscopy, (b) Magnified view of the spectroscopy around 2.88 GHz. These 
pictures clearly show that every peak has a finite energy shift if we add the detuning. 
Here, we use the same parameters as in Fig. [TT] 
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1 g g 2 S J 2 5 

+ 7! + 4 (s 2 + J2 ) 3/2 (s 2 + J2 ) 3/2 ' T 

1 J r, A •/-■> . . , 

+ v / 2 V »^ 5 4 ( 9 2 + J2 ) 3/2 (» 2 + ^ 2 ) 3/2 ' W 

These results show that, as the detuning increases, the energy difference between the 
excited states and the ground state would increase. The amount of change depends 
on g, J, and 5 (shown in Fig JT5]l . It is worth mentioning that the middle peaks has a 
smaller energy shift than that of the other two side peaks for g J. 

From the spectroscopy calculated by MHOM, we have also confirmed the fact that 
the frequency of middle peak can be shifted by the detuning (See Figjl5]and Fig JlBl) . 
Also, the left and right peaks have an energy shift for the detuning, and the energy shift 
is much larger than that of the middle peak. ( See the FigH5]and Fig JTfil) 
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